Time delay in QSO 0957+561 from 1984-99 optical data 
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ABSTRACT 

Photometric optical data of QSO 0957+561 covering the period 1984-99 are analyzed 
to discern between the two values of the time delay (417 and 424 days) mostly accepted 
in the recent literature. The observations, performed by groups from three different 
institutions — Princeton University, Harvard-Smithsonian Center for Astrophysics, and 
Instituto de Astroffsica de Canarias — and including new unpublished 1998-9 data from 
the IAC80 Telescope, were obtained in five filters {V, R, I, g, and r). The different 
light curves have been divided into observational seasons and two restriction have been 
applied to better calculate the time delay: (i) points with a strange photometric behavior 
have been removed; and (ii) data sets without large gaps have been selected. Simulated 
data were generated to test several numerical methods intended to compute the time 
delay (Atab)- The methods giving the best results — the discrete correlation function, 
(5-square, z-transformed discrete correlation function, and linear interpolation — were 
then applied to real data. A first analysis of the 23 different time delays derived from 
each technique shows that Atab must be into the interval 420-424 days. From our 
statistical study, a most probable value of Atab = 422.6 ± 0.6 days is inferred. 

Subject headings: quasars: individual: Q0957+561 — cosmology: gravitational lensing — 
photometry: DAOPHOT — methods: data analysis 



1. Introduction 

The first discovered gravitational lens system, QSO 0957+561 (Walsh, Carswell, & Weymann 
1979), has been the subject of a continuous and exhaustive monitoring in several bands since 1979. 
The special characteristics of this system made it very attractive for time delay determinations, and 
different values for this quantity were presented during the 1980s: At = 566 ± 37 days (Florentin- 
Nielsen 1984); At = 376 ± 37 days (Schild & Cholfin 1986); At = 657 ± 73 days (Gondhalekar et 
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al. 1986); At = 478 ± 73 days (Lehar, Hewitt, &; Roberts 1989). As can be seen, there was wide 

dispersion in the results obtained by different groups. However, the monitoring campaigns carried 
out during the early 1990s led to quite an odd situation, as all the results concentrate around two 
different values for the time delay: Ar ~ 420 days and At ~ 510 days. Calculations leading to 
the first value were presented by Vanderriest et al. (1989, Ar = 415 it 20 days); Schild (1990, 
At = 404 ± 10 days); and Pelt et al. (1994, At = 415 ± 32 days in the B band and At = 409 ± 23 
days in radio). On the other hand, a value close to 510 days was obtained by Beskin & Oknyanskij 
(1992, Ar = 522 it 15 days in the B band and Ar = 515 it 15 days in the r band); Roberts et al. 
(1991, Ar = 515 ± 37 days); and Press, Rybicki, & Hewitt (1992a, Ar = 536l}2 days in the B 
band; 1992b, Ar = 540 ± 12 days in B+radio). 

This situation abruptly changed when Kundic et al. (1995) presented their observations in the 
g and r bands. A sharp drop appeared in 1994 December and could be used to discern between the 
"long" (510 days) and the "short" (420 days) time delay, provided continuous monitoring of QSO 
0957+561 was carried out in the first half of 1996. This monitoring was performed, and the long 
time delay was rejected (Oscoz et al. 1996; Kundic et al. 1997). The controversy regarding the 
time delay seemed to be finally solved. However, the results appearing in the literature since 1995 
concentrate again around two values, 417 and 424 days. These results are summarized as follows: 
Ar = 423 ± 6 days (Pelt et al. 1996); Ar = 417 ± 3 days (Kundic et al. 1997); Ar = 424 ± 3 
days (Oscoz et al. 1997); Ar = 425 ± 17 (Pijpers 1997); Ar = 416.3 ± 1.7 days (Pelt et al. 1998); 
At = 425 ± 4 days (Serra-Ricart et al. 1999); Ar = 417.4 days (Colley & Schild 2000). 

This difference is irrelevant in the Hubble constant calculations, as the uncertainty introduced 
by the time delay is much lower than the uncertainty given by other factors (e.g., the main lens 
galaxy's velocity dispersion or the lens modeling). However, the most accurate time delay should 
be used in the search for possible very rapid microlensing events in QSO 0957+561, and a week's 
difference in Ar could lead to the detection of false events or failure to detect real ones. 

In this paper we have compiled photometric data of QSO 0957+561 from three different ob- 
serving groups covering the period 1984-99 to obtain an estimate of the time delay by means of 
several statistical methods. This includes new unpublished data corresponding to the last obser- 
vational campaign (1998-9) at the IAC80 Telescope. The data sets are presented in §2, while the 
methods for obtaining the time delay appear in §3. A first check of the goodness of these methods 
applied to simulated data is calculated in §4, and the best techniques are applied in §5 to real data. 
Finally, a discussion of our results appears in Section 6. 



2. Selected data sets 

Several monitoring campaigns of QSO 0957+561 in different bands have been performed since 
1979. However, to obtain the time delay, only the observations obtained by groups from three differ- 
ent institutions will be considered here: Princeton University (hereafter PU), Harvard-Smithsonian 
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Center for Astrophysics (hereafter CfA), and the Instituto de Astrofi'sica de Canarias (hereafter 
lAC). 

2.1. PU data 

Images were obtained at the Apache Point Observatory 3.5 m telescope in the g and r bands. 
QSO 0957+561 was monitored during several observational campaigns, although only data corre- 
sponding to the first two seasons have been published (Kundic et al. 1995, 1997): (i) from 1994 
December to 1995 May; and (ii) from late 1995 to mid 1996. The resulting data comprise 51+46 
g-hand points and 54+46 r-band points. The light curves were calculated via aperture photometry, 
and have neither large error bars nor significant gaps. Their main characteristic is the presence of a 
sharp drop of about 0.1 mag in the A-component in late 1994 December, very useful for time-delay 
calculations. 

2.2. CfA data 

This data set is the largest ever obtained for a gravitational lens system. It consists of 1069 
brightness measurements in the R band, from late 1979 to mid 1996. The observations correspond- 
ing to the period 1979-89 were obtained at the Whipple Observatory 0.61 m telescope, while the 
remaining data were obtained with a 1.2 m telescope (Schild &; Thomson 1995, and references 
therein). The reduction procedure followed a basic aperture photometry scheme (although a new 
automated photometry reduction code is now being applied by the authors, the results do not 
substantially differ from the "old" photometry). The error bars arc not large, with the exception 
of the first five years. The main problem with this data set is the scarcity of observations during 
the first 1800 days (81 brightness measurements). Moreover, those points coincide with the largest 
error bars. So, we will consider the data from mid 1984 for time delay calculations. 

2.3. lAC data 

Lens monitoring was performed in four consecutive seasons (1996 February to June, 1996 
October to 1997 July, 1997 October to 1998 July, and 1998 October to 1999 July) using the CCD 
camera of the 82 cm IAC-80 Telescope (IAC80 hereafter), sited at the lAC's Teide Observatory 
(Tenerife, Canary Islands, Spain). A Thomson 1024x1024 chip was used, offering a field of nearly 
T.b. Standard VRI broad-band filters were used for the observations, corresponding fairly closely 
to the Landolt system (Landolt 1992). The final data set comprises 172 point in the V band, 
301 points in the R band, and 112 points in the / band. Accurate photometry was obtained by 
simultaneously fitting a stellar two-dimensional profile on each component by means of DAOPHOT 
software (see details in Serra-Ricart et al. 1999). A new, completely automatic IRAF task has been 
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developed demonstrating, using a sample of simulated data, that the proposed method can achieve 
high-precision photometry. However, the errors bars obtained for IAC80 data are slightly larger 
than those of PU and CfA data. This could be explained by a decrease in chip sensibility due to 
the age of the CCD. In order to assess the reliability of our method using real data, simultaneous 
observations of QSO 0957+561 were undertaken on 1999 February 19 by using the IAC80 and 
the 2.5 m Nordic Optical Telescope (NOT hereafter) sited at the lAC's Roque de los Muchachos 
Observatory (La Palma, Canary Islands, Spain). The final reduced results arc presented in Figure 
1 and Table 1 (photometric errors for comparison star H and D arc also included). The NOT 
light-curve errors (a few millimagnitudes) are much lower than the IAC80 ones, and this difference 
could be explained in terms of the following: i) the NOT has a larger aperture than the IAC80, 
and ii) the NOT also has a better CCD chip. However, the good agreement between the two curves 
demonstrates that our reduction method works with high degree of accuracy. The photometric 
data are available at URL http://www.iac.es/lent. 



The large amount of data described in §2 adds new biases (different telescopes, filters, reduction 
processes, and behavior) to the inherent difficulty to analyse discrete, unevenly sampled temporary 
series. These facts leaded us to employ several statistical methods to calculate the time delay 
to increase the robustness of the results thus obtained. As a first step, several techniques will 
be checked by using simulated data: the discrete correlation function, dispersion spectra, 5"^, 5"^ 
modified, linear interpolation, and the z-transformed discrete correlation function. 



The DCF (Edelson & Krolik 1988) is a technique valid for any physical quantity that is observed 
to vary in time. For two discrete data trains, Ai and Bj, the formula representing their DCF is 



averaging over the M pairs for which r — a < A.tij < t + a, a, e/j, and a/ being the bin semi-size, 
the measurement error associated with the data set fc, and the standard deviation, respectively. 
The maximum of the DCF is identified with the time delay. 



The data model (Pelt et al. 1996) consists of two time series, = q{ti) + eAiU) , i = 1,. . . ,Na, 
and Bj = q(tj — ATBA) + l{tj) + ^B{tj), j = 1, - ■ ■ ,Nb, where q{t) represents the intrinsic variability 



3. Different methods to obtain time delay estimates 



3.1. Discrete correlation function (DCF) 




(1) 



3.2. Dispersion spectra 
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of the quasar, l{t) accounts for the difference in magnitudes plus additional variability in time due 

to microlensing, and e^(t) and es^t) are observational errors. This two series are combined into 
one, C, for every fixed combination [r, l{t)] by taking all values of A as they are and correcting the 
B data by l{t) and shifting them by r. The dispersion spectra that will be used here are represented 

by 



^'^ i{t) S^''^ W G 



^4,it = mm -7TX— , (2) 



where the Wn,m are statistical weights of the combined light curves, and Gn,m = 1 when C„ and Cm 
come from different curves {A or B) and otherwise. From eq. (2), we can consider two different 
approximations depending on the definition of Srtjn'- (i) Snjn = 1 if \tn — tm| < <^ and otherwise; 

(2) ' ' 

and (ii) Sn^m = 1 — \tn — tm\/5 if \tn — Un\ < ^ and otherwise, 5 being the maximum distance 
between two observations which can be considered as nearby. The minimum value of eq. (2) is 
assumed as the time delay. 



3.3. d-square 

The (5^ method (Serra-Ricart et al. 1999) makes use of the similarity between the discrete 
autocorrelation function (DAC) of the light curve of one of the components and the A-B discrete 
correlation function (DCF). Prom the DAC and DCF functions, one can define a function 

/ 1 \ ^ 

^-(^) = U E [I^CF(r,) - DAC(r, - e)f (3) 

^ ^ i=i 

for every fixed value 6 (days), where Si = 1 when both the DCF and DAC are defined at Tj and 
Ti — 9, respectively, and otherwise. The most probable value for the time delay should correspond 
to the minimum of this function. 



3.4. modified 

This modification of the 6^ method was suggested by Schild (1999, private communication). 
It consists in comparing the DAC and DCF curves by taking their ratio instead of by calculating 
their difference. So, the final equation to obtain the time delay is 

«('')=(^)i:«- 

^ ^ i=i 



DCF(r,) 
DAC(rj -e) 



(4) 
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3.5. Linear interpolation (LI) 

The linear method is similar to that suggested by Kundic et al. (1997). One of the two light 
curves (hereafter light curve 1) is selected, and the linear interpolation of data and their errors is 
considered as reference. The other light curve (hereafter light curve 2) is then shifted in magnitude 
by jTist the difference between the means of both light curves. After that, light curve 2 is shifted in 
time and chi-square per degree of freedom (x^) for each time delay is calculated. The number of 
degrees of freedom is equal to the number of points of the light curve 2 in the overlapping interval 
minus 2 (because we are fitting shifts in magnitude and time). The time that minimizes is taken 
as a provisional time delay. 

This procedure is followed by using as reference both the A- and U-component light curves, 
selecting then the time delay closest to 421 days (the intermediate point between 417 and 425 days). 
The uneven sampling of the light curves usually leads to a better time delay taking as a reference 
one of the two light curves. 

3.6. ^-transformed discrete correlation function (ZDCF) 

The ZDCF (Alexander 1997) is a new method for estimating the cross-correlation function 
(CCF) of sparse, unevenly sampled light curves. Fisher's z-transform of the linear correlation 
coefficient, r, is used to estimate the confidence level of a measured correlation. This technique 
attempts to correct the biases that affect the original DCF by using equal-population binning. The 
ZDCF involves three steps: 

(i) All possible pairs of observations, {ai,bj}, are sorted according to their time-lag, ti — tj, 
and binned into equal population bins of at least 11 pairs. Multiple occurrences of the same point 
in a bin are discarded so that each point appears only once per bin. 

(ii) Each bin is assigned its mean time-lag and the intervals above and below the mean that 
contain la of the points each. 

(iii) The correlation coefficients of the bins are calculated and ^-transformed. The error is 
calculated in 2;-space and transformed back to r-space. 

The time-lag corresponding to the maximum value of the ZDF is assumed as the time delay 
between both components. 

4. Analysis of the different methods by using simulated light curves 

The application of the statistical methods described in §3 to simulated data sets can serve to 
check the validity of their results under different conditions, always with discrete and irregularly 
sampled data sets. The six selected data sets are quite similar to those presented in Serra-Ricart et 
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al. (1999), where several sets of simulated photometric data with similar irregularity in the obser- 
vations (time distribution of the data), magnitudes, and error bars (i.e.. as large as or even larger 
than those of PU, CfA, and lAC; the worst situation is selected) to that of the lAC observations 
were created. 

First, a set of dates, .x,, between 1800 and 2000 (TJD = JD-2449000) approximately, was 
generated with a pseudo-random separation, taken from a uniform distribution between and 5 
days. These data were alternatively separated in two time series, corresponding to A- and B- 
component light curves, this last curve being shifted by 420 days to simulate the existence of a time 
delay. A first magnitude was calculated for each date with the relationship yi = F{xi) [see below for 
the different shapes of F(xj)]. The probability of measuring a value y for each Xj is proportional to 
Qi-iy-Vi) /2o-J^ hence characterized by ctj, or, equivalently, the variable d = y — yi is distributed 
as e^~'^^ ^"^"i^ . A ai taking pseudo-random values between 0.01 and 0.03 was generated for each Xi. 
Prom here the quantities di, pseudo-random numbers obtained from a normal Gaussian distribution 
with zero mean and standard deviation, o"j, were calculated, allowing them to adopt positive or 
negative values. Pinally, the magnitude was generated from the equation yg = F{xi) + di = yi + di, 
with an error bar of crj. The A component was made brighter by adding 0.1 to the magnitudes of 
the B component. 

The first selected function was: 

Fl: y = 17.17 + 0.5 e-"-^-^ sin(/), where / = ■ (5) 

This function represents light curves in which a sharp event similar to that reported by Kundic et 
al. (1995) appears. 

The second function is 

P2 : y = 17.2 + 0.1 sin(/) sin(4/), where / = ^ • (6) 

In this case, the light curves present several maxima and minima, although none of them is clearly 
remarkable. 

An additional function, consistent with the actual variability of Q0957+561, was created (the 
lAC observational data from 97-98 seasons, were selected as reference). The light curves were then 
fitted by the function 

F3:y = 17.07 - O.lGe^, where / = 'y^^^^ (7) 

m being the mean of the TJD in the selected range and s its standard deviation. The resulting 
simulated data show a lower variability to that obtained from Fl and F2. 

Besides the comparison between the A- and i?-component for the three data sets, an additional 
test was performed by removing some data of the ^-component. By doing this, we want to simulate 
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the data sets corresponding to certain periods in which no vahd data could be obtained during 
several days (bad weather, problems with the telescope, etc.). 

The time delay corresponding to the different statistical techniques were first calculated by 
allowing their free parameters to vary. The best results were obtained with a a = 5 days for the 
DCF, a = 2 days for D| ^ and 2, and 5 = 10 days for 5^ and 6"^^. In the DCF, d'^ and 5^ cases, 
the results were quite similar for the parameters varying between 4 and 12 days. However, slightly 
better values were obtained for a = 5 days in the DCF case, and for ^ = 10 days in the S'^ and 
5^ cases, so these quantities were selected. The objective in this paper is to analyze the different 
methods with several data sets by using the same conditions. An analysis of the best method 
and/or parameters to be used with a particular data set will be presented in a future paper. The 
uncertainties in the time-delay estimates were computed by generating 1000 bootstrap samples and 
applying the statistical methods in each case. 

Some interesting consequences can be derived from the results in Table 2. As expected, the 
best estimate of the time delay is always obtained for the function Fl, i.e., when a sharp drop 
similar to the one appearing in §2.1 is present in the light curves. All the methods give a good 
value for the delay. On the other hand, in the F2 and F3 cases the error bar are too large for 
S^, ^4,1, and D4 2 when compared with the results obtained with the other methods. So only the 
DCF, (5^, ZDCF, and LI techniques will be applied to real data calculations, similar in most of the 
cases to the functions F2 and F3. Notice that the true time delay is always within the error bar, 
even when a large gap is present in the light curves. 

5. Time delay from real data 

Prior to applying the different methods to calculate the time delay from real data, some 
considerations have to be taken into account. First, data should be checked to eliminate incon- 
sistent measurements. This modification of the raw data, based on a suggestion by Falco (1997, 
private communication), takes account of the possible existence of strong and simultaneous (not 
time-shifted) variations of some data point in both components. The inclusion of such "strange" 
brightness records in the final data sets probably originated from failures in the CCD or from bad 
weather conditions, creates artificial peaks or valleys in the light curve of one of the components. 
These maxima/minima have no importance when a sharp change in the behavior of the quasar 
is being analyzed, but can lead to a completely wrong time delay estimate when dealing with a 
smoother season. To avoid these abnormal observations we have removed the points with a simul- 
taneous difference in magnitude in both components as compared with the previous and following 
records larger than 2.5 times their error bar. This was done by considering only those points with 
a difference in the observation dates of less than 10 days. The resulting data sets will be named 
bad point free (BPF). The 301 points of the lAC data set in the R band are reduced to 289 with 
this restriction, while less than 60 of the 1069 CfA observations have to be removed. Finally, no 
brightness measurement of the PU data seems to be wrong. However, the BPF restriction could be 
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applied neither to the lAC /- and F-band data nor to the CfA 89-90 and 91-92 i2-band data due 
to their low number of points, with a minimum distance between neighboring points more than 10 
days in most of the cases. 

Another drawback when dealing with real data is the impossibility of observing the system 
during certain months of the year and the consequent lack of suitable edges. Once it is stated, see 
§1, that the rough value of the time delay is around 420 days, the comparison between the A and 
B components should be made by previously selecting a "clean" data set (CD, see Serra-Ricart et 
al. 1999), i.e., homogeneous monitoring of both images during two active and clear (free from large 
gaps) epochs separated by ~ 420 days. 

Finally, both types of corrections, BPF and CD, were combined to obtain the definitive data 
sets used in this paper. To see the importance of applying the CD-BPF corrections, we mention 
two extreme examples: 1) At = 398 it 11 days with the raw lAC data corresponding to the 96-7 
season, while At = 428 it 9 days with the CD-BPF approximation; and (2) At = 384 it 5 days 
with the original 93-4 CfA data and At = 423 zt 2 days with the CD-BPF approximation. 

To summarize, Monte Carlo calculations were applied to the four techniques (DCF, (5^, ZDCF, 
and LI) with the CD and the CD-BPF restrictions. The PU, lAC, and CfA data sets were divided 
into observational seasons, leading to the 23 different time-delay estimates per method appearing 
in Tables 3 to 6. The results obtained from the PU data in both the g and r filter are presented in 
Table 3. Notice that two of the methods (DCF and LI) employed here are also used by Kundic et 
al. (1997). Our results are quite similar to those obtained by these authors, and, moreover, they 
are always into their error. The small differences come from the selection of clean data sets. Tables 
4 and 5 offer the time delay for the lAC data in the R, and in the I and V bands, respectively. 
Finally, Table 6 was obtained from the CfA R band data. The analysis of these delays are done by 
considering the CD-BPF quantities, except in those cases in which only the CD results could be 
obtained. The quantities obtained for the CfA 92-3 season with the DCF and methods can be 
discarded, as they appear to be clearly inconsistent (394 zt 1 and 403 ± 1 days, respectively). 

6. Discussion and conclusions 

A first step in discussing the value of Atab consists in computing, for each of the techniques 
(DCF, (5^, ZDCF, and LI), the number of occurrences of each value of the time delay. These 
quantities, obtained from Tables 3 to 6, are depicted by the black lines in Figure 2: DCF (top 
panel, left); 5^ (top panel, right); ZDCF (bottom panel, left); and LI (bottom panel, right). (The 
different values of Atab have been grouped into two-day bins). As can be seen, two remarkable 
characteristics can be deduced from Figure 2: (i) there is a small dispersion in the delays, as most 
of them are in the interval 415-428 days; and (ii) the centroids of the histograms, given by the 
average of the time delays derived in Tables 3 to 6, are always in the interval 420-424 days. These 
last quantities are represented in Figure 2 (open circles) together with their uncertainty (r.m.s., see 
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Table 7 and discussion below). Note that the largest peak of each histogram coincides with this 
average value, except for the ZDCF technique. In the histogram corresponding to 5^, the maximum 
corresponds to the values 421-424, while the average is 421.8 it 1.3 days. The DCF panel docs not 
lead to a clear time delay, with two maxima in the intervals 417-420 and 423-424 days. The average 
here is given by 423.3 it 1.4 days. In the case of LI, the peak is placed at At^b = 423-424 days, 
and the average is 424.3 it 1.2 days. Finally, the ZDCF method has a maximum around 423-424 
days. However, the average is 420.6 ib 1.1 days, which is slightly different. The red histograms in 
Figure 2 represent the total number of occurrences obtained by adding the results from the four 
techniques. Its center is again over 420 days, giving a maximum of 423-424 days. 

To complement these calculations, which have been done by fixing the method and computing 
the probability of appearance of each delay, we can now represent the number of times each value 
appears for each data set, independently of the method employed. Four different data sets have 
been selected: (i) PU r and g filters, Figure 3a; (ii) lAC R data. Figure 3b; (in) lAC V and I 
records. Figure 3c; and (iv) CfA R data. Figure 3d. In this case, the centroid of the distributions, 
represented by open circles in Figure 3, is again over 420 days: 421.4 it 1.1 days, 423.7 ± 1.3 days, 
423.7 ± 1.2 days, and 421.8 ± 1.0, for PU r and g, lAC R, lAC / and V, and CfA R, respectively. 
A first positive consequence of the PU results is their extremely short dispersion indicating the 
goodness of the data and the presence of the sharp event. The maximum here is placed between 
419 and 424 days. However, the clearest peak appears from the lAC R values around 423-424 days. 
The lAC / and V panel shows more dispersion, probably due to the higher error bars of the light 
curves. There is not a unique maximum here, as two peaks appear around 423-424 and 427-428 
days. The highest dispersion in the results can be seen in Figure 3d, corresponding to the CfA 
data. The maximum would in this case be in the interval 417-422 days. A remarkable result here 
is that the average coincides with the maxima in all the panels. Once again, the red histogram 
gives the total number of occurrences. 

The combination of the results derived from Figures 2 and 3 (centroids and maxima of the 
distributions) support a 1\tab in the range of 420-424 days, although a time delay of around 
417 days cannot be totally discarded. It is important to remark that the results are the same 
independently of the technique employed or of the data set selected. 

The time delay between A and B components of Q0957+561 does not depend on the filter 
(as it is an acromatic effect) and/or the time (different campaigns), so it should be possible to 
merge the different sample results. A very important point is to assess the statistical reliability 
of the different delay calculation methods in order to estimate final delay errors. Several statistics 
were calculated. Mean Delay {MD = J2^=i ArABi / N) , Mean Error {ME = Y^^^^ei/N, with 
individual errors ), and Dispersion {DI = [J2Zii^^ABi - MDf/{N - 1)]^/^). If the error estimate 
is correct, then ME DI, and the final error for the time delay will be given by the r.m.s. 
{[T.f=i{^^ABi - MDf/{N *{N- l))]^/2) (see Eadie et al. 1971). Tables 7 to 10 (see below) show 
the final results for the four methods. In all cases, within the statistical errors, good agreement is 
found between the mean error and the dispersion. 
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When this procedure is apphed to the different time delays given by each technique, one 
obtains the results in the first four rows of Table 7, where the mean values of the time delay 
and the uncertainties (r.m.s.) are shown. According to these calculations, the definitive time 
delay would be in the interval 420(ZDCF)-424(LI) days, with an uncertainty below 1.4 days. This 
interval coincides with that derived from the analysis of Figures 2 and 3. When this calculation is 
done with only the R band results (Table 8), the results are almost identical. 

The different analyses performed until now have been done considering all the results obtained 
in Tables 3 to 6. However, the error bars of some of these values exclude the interval 420-424 
days, where, as shown before, there is the highest probability of finding the right Atab- The mean 
delays and uncertainties obtained when these values are removed (first four rows of Table 9) are 
very similar to those of Table 7. Atab is now restricted to the interval 420.8-423.2 days, with an 
uncertainty below 1.3 days. Once again, the results derived from the i?-filter data are almost the 
same (Table 10). 

The validity of these statistical results leaded us to repeat the calculations of Tables 7 to 10, 
but this time having into account all the delays, i.e., without considering the method (DCF, S'^, 
ZDCF, or LI). This allows us to have four different time delays per year in most of the occasions, 
and so, a larger amount of data for the statistical analysis. The results appear in the last row 
of Tables 7 to 10, although only the values of Tables 7 and 9, Atab = 422.6 ± 0.6 days and 
Atab = 422.0 it 0.6 days, respectively, will be considered, as they were obtained from a larger 
amount of data. Adopting a conservative point of view, we will select the quantity with the higher 
uncertainty, Atab = 422.6 it 0.6 days, as the final time delay. 

Different treatments of the time delays obtained in §5 have been performed. The analysis of 
these results always points to a time delay in the interval 420-424 days. None of these methods 
clearly favors the values 416-418 days as the right time delay. Moreover, the averages of the 
quantities of Tables 7 and 9 give values around 422 days, coinciding with the maxima obtained 
in Figures 2 and 3. Assuming this value as the time delay between the A and B components of 
Q0957+561, let us check which of the results given in Tables 3 to 6 include 422 days in their the 
error bars. Figure 4 offers the number of times, written as a percentage of the total number of time 
delays obtained for each method, that each value of Atab is included in these error bars. This 
probability is the same, 86%, for DCF, and ZDCF, and 74% for LI. In this case, the DCF and (5^ 
methods show the highest probability for a Atab of 421-422 days, while the peak in the LI curve is 
placed at 423 days. A wider maximum is obtained in the case of ZDCF, with the same probability 
between 418 and 421 days. 

As can be seen, not all the seasons in the different data sets are fully appropriate for calculating 
the time delay. However, our intention was to analyze all the data in the three different data sets 
with four different methods, and finally to restrict the uncertainty in the calculation of the time 
delay. The statistical treatment of all the results confirms a time delay of 422 days. 
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Fig. 1. — Light curves for A and B images of QSO 0957+561 in the R band obtained on the night 
of 1999 February 19. IAC80 data are represented by boxes and NOT data are indicated by circles. 
One-sigma error bars are indicated. See text for more details. 

Fig. 2. — Number of times (black lines) that each value of the time delay appears, obtained from 
the results of Tables 3 to 6. The four different techniques have been considered: a) DCF, b) (5^, c) 
ZDCF, d) LI. The open circles represent the average value of the delays for each method, and the 
red histograms represent the sum of the values given by each method. 

Fig. 3. — Number of times (black lines) that each value of the time delay appears for each data set 
from the values given in Tables 3 to 6. Four data sets are represented: a) PU r and g, b) lAC R, 
c) lAC / and V, and d) CfA R. The open circles represent the average value of the delays for each 
data set, and the red histograms represent the total sum of the values. 

Fig. 4. — Percentage of times that a value of the time delay is included in the error bars of the 
results given in Tables 3 to 6: black = DCF, red = 6'^, green = ZDCF, and blue = LI. 
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Fig. 4 
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Table 1. R photometric errors 



Light curve^ 


NOT*^ 


IAC80 


A 


0.005 


0.016 


B 


0.005 


0.014 


H-D 


0.002 


0.005 



^See section 1 for details. 

''Error values are expressed 
in magnitudes. 
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Table 2: "Time delays" obtained from the application of a Monte Carlo algorithm with the different 
techniques to the six simulated data sets (see main text for details). The "true" time delay was 
420 days. 



Method 




Fl 


Fl (gap in A) 


F2 


F2 (gap in A) 


F3 


F3 (gap in A) 


DCF (q 


= 5) 


421±2 


421±4 


419±1 


420±2 


415±9 


412±9 


{6 = 


10) 


422±2 


421±2 


419±1 


420±1 


420±5 


417±5 


ZDCF 




420±3 


420±3 


421±1 


420±7 


418±14 


416±13 


LI 




421±2 


421±2 


420±1 


421±1 


419±12 


418±13 


Si iS = 


10) 


422±9 


422±13 


425±30 


416±34 


403±34 


413±29 


Dli (a 


= 2) 


421 ±3 


421±4 


421±12 


416±20 


422±27 


414±25 


Dh (« 


= 2) 


421 ±2 


421±3 


423±15 


415±26 


420±26 


412±21 
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Table 3. Time delays obtained from the application of a Monte Carlo algorithm with four 

different methods to the PU r and g data 



Filter Years 


Method 


Points 


Delay 






(CD) 


(CD) 




DCF 




422±1 


g A (94-5); B (95-6) 


52 


42; 39 


423±1 




ZDCF 




416±2 




LI 




419±5 




DCF 




426±5 


r A (94-5); B (95-6) 


52 


41; 41 


423±2 




ZDCF 




420±8 




LI 




422±3 



-23- 



Table 4. Time delays obtained from the application of a Monte Carlo algorithm with four 

different methods to the lAC R data 



Years 


Method 


Points 


Delay 


Points 


Delay 






(CD) 


(CD) 


(CD-BPF) 


(CD-BPF) 




DCF 






428±12 




425±11 


A (96-8); B (97-9) 




182; 


220 


421±5 


173; 212 


421±4 




ZDCF 






418±13 




426±13 




LI 






421±16 




424±17 




DCF 






432±8 




436±8 


A (96); B (96-7) 


6^ 


31; 


36 


425±8 


28; 33 


428ib9 




ZDCF 






422±13 




424±14 




LI 






424±7 




424±7 




DCF 






436±10 




424±13 


A (96-7); B (97-8) 




46; 


86 


431±4 


45; 84 


425±4 




ZDCF 






429±13 




424±12 




LI 






415±16 




426±16 




DCF 






414±12 




416±12 


A (97-8); B (98-9) 




79; 


62 


414±10 


78; 61 


413±8 




ZDCF 






420±12 




420±13 




LI 






427±13 




423±13 
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Table 5. Time delays obtained from the application of a Monte Carlo algorithm with four 



different methods to the lAC / and V data 


Filter 


Years 


Method 


Points 
(CD) 


Delay 
(CD) 


I 


A (96-8); B (97-9) 


DCF 

6^ 
ZDCF 

LI 


65; 87 


415±13 
423±12 

417±14 
419±9 


I 


A (96); B (96-7) 


DCF 

ZDCF 
LI 


19; 18 


424±16 
424±16 

430±3 


I 


A (96-7); B (97-8) 


DCF 

5^ 
ZDCF 

LI 


18; 31 


419±9 
422±7 

427±7 


I 


A (97-8); B (98-9) 


DCF 

6^ 
ZDCF 

LI 


28; 39 


414±15 
417±18 

421±14 
427±11 


V 


A (96-8); B (97-9) 


DCF 

6^ 
ZDCF 

LI 


68; 149 


432±11 
429±11 

418±9 
428±14 


V 


A (97-8); B (98-9) 


DCF 

52 
ZDCF 
LI 


29; 84 


433±7 

432±6 
428±9 
423±3 
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Table 6. Time delays obtained from the application of a Monte Carlo algorithm with four 



different methods to 


the CfA R data 




Years 


Method 


Points 
(CD) 


Delay 
(CD) 


Points 
(CD-BPF) 


Delay 
(CD-BPF) 


A (84-5); B (85-6) 


DCF 

52 
ZDCF 

LI 


40; 49 


420±6 
420±5 
429±7 
432±5 


39; 47 


419±4 
421±6 
428±7 
431±6 


A (85-6); B (86-7) 


DCF 

52 
ZDCF 

LI 


49; 74 


419±10 
410±13 
420±9 
409±14 


46; 72 


417±12 
407±15 
417±12 
431±8 


A (86-7); B (87-8) 


DCF 

52 
ZDCF 

LI 


53; 58 


431±12 
434±10 
429±15 
424±9 


52; 58 


428±13 

418±17 
430±12 
420±5 


A (87-8); B (88-9) 


DCF 

52 
ZDCF 

LI 


60; 53 


425±12 
423±17 
419±17 
421±11 


58; 53 


428±13 
422±18 
411±13 
422±12 


A (88-9); B (89-90) 


DCF 

52 
ZDCF 

LI 


23; 19 


417±21 
420itl0 
412±8 
421±2 






A (89-90); B (90-1) 


DCF 

52 
ZDCF 

LI 


40; 38 


440±5 
407±7 
417±6 
410±6 


40; 34 


424±19 
420±3 
420±7 
411±3 


A (90-1); B (91-2) 


DCF 

52 
ZDCF 

LI 


15; 29 


435±8 
433±18 

423±19 
426±6 






A (91-2); B (92-3) 


DCF 

52 
ZDCF 

LI 


14; 72 


394±0 
402±1 
423±12 
436±9 


14; 67 


394±1 
403itl 
416±13 
437±8 


A (92-3); B (93-4) 


DCF 

52 
ZDCF 

LI 


70; 98 


411±5 

423±2 
411±11 
423±13 


63; 93 


419±6 

423±2 
424±15 
433±3 


A (93-4); B (94-5) 


DCF 

52 
ZDCF 

LI 


83; 111 


422±4 
421±2 

422±6 
421±4 


78; 95 


422±4 
421±2 

422±7 
418±6 


A (94-5); B (95-6) 


DCF 

52 
ZDCF 

LI 


101; 66 


415±6 
411±2 
416±2 
421±12 


87; 57 


418±7 
415±3 
416±5 
418±5 
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Table 7. Statistical results obtained from the 23 time delays calculated for each method 



Method 


Mean Delay 


Mean Error 


Dispersion 


r.m.s. 


DCF 


423.3 


10 


7 


1.4 




421.8 


9 


6 


1.3 


ZDCF 


420.6 


11 


5 


1.1 


LI 


424.3 


7 


6 


1.2 


Total 


422.6 


9 


6 


0.6 



Table 8. As in Table 7 but only considering R filter time delays 



Method 


Mean Delay 


Mean Error 


Dispersion 


r.m.s. 


DCF 


423.4 


11 


6 


1.7 




420.5 


8 


6 


1.7 


ZDCF 


420.9 


11 


6 


1.4 


LI 


424.3 


8 


7 


1.7 


Total 


422.3 


10 


6 


0.8 
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Table 9. Statistical results obtained only considering time delays and error bars than include 
part or the whole interval given by the results of Tabic 6 (420-424 days) 



Method 


Mean Delay 


Mean Error 


Dispersion 


r.m.s. 


DCF 


422.2 


11 


6 


1.3 


5^ 


421.6 


9 


6 


1.2 


ZDCF 


420.8 


11 


5 


1.2 


LI 


423.2 


8 


4 


0.9 


Total 


422.0 


10 


5 


0.6 



Table 10. As in Table 9 but only considering R filter time delays 



Method 


Mean Delay 


Mean Error 


Dispersion 


r.m.s. 


DCF 


422.5 


11 


6 


1.6 




420.9 


9 


6 


1.8 


ZDCF 


420.9 


11 


6 


1.4 


LI 


423.0 


8 


4 


1.2 


Total 


421.7 


10 


5 


0.8 



